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Abstract 

An  interacting  lattice  gas  model  is  used  to  study  flow  of  immiscible  components  A  and  B  (molecular  weights  Ma  and 
Mb,Ma<^b)  by  Monte  Carlo  simulations.  Concentration  gradients  and  hydrostatic  pressure  bias  (H)  drive  these 
constituents  from  their  source  at  the  bottom  against  gravitational  sedimentation  in  an  effective  medium.  Response  of  their 
flux  densities  (/^  to  the  hydrostatic  bias  H  are  examined.  If  both  constituents  are  released  with  equal  probabilities  (a 
non-interacting  source),  their  flux  densities  respond  linearly  to  bias  with  jA>jB  except  at  the  extreme  bias  H  — ►  1  where 
Ja  Jb'  f^^ow  response  becomes  complex  if  the  constituents  from  their  source  are  released  according  to  their  current 
lattice  concentrations  (an  interacting  source):  a  crossover  occurs  from at  low  bias  (H  <  0.4)  to  jB>jA  higher  bias 
(//>0.4).  Constituent  with  the  lower  molecular  weight  {A)  responds  linearly  on  increasing  the  bias  except  at  very  high  bias 
(H  >  0.8)  where  the  response  becomes  negative.  The  heavier  component  {B)  responds  non-linearly:  a  high  response  at  low 
values  of  H  is  followed  by  a  linear  response  before  the  onset  of  eruptive  response  at  high  range  of  H.  The  volatility 
parameter  diverges  as  eruption  occurs  at  //  -►  1 . 

Published  by  Elsevier  B.V. 
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1.  Introduction 

Flow  patterns  and  morphological  evolution  in  the  sub-bottom  region  (at  and  below)  of  the  ocean  floor  are 
complex  due  to  dynamics  of  interacting  components  with  diverse  characteristics  such  as  granular  matters 
(sand,  silt,  etc.),  mixtures  of  fluids,  and  gas  [1-7].  The  response  of  these  constituents  to  temperature,  pressure, 
and  concentration  gradients  lead  to  a  variety  of  correlated,  local  and  global  phenomena  involving  flow, 
extrusion,  evolution,  and  settling  of  constituents  with  a  range  of  relaxation  time  and  length  scales  [4,8-10]. 
Understanding  of  structure  and  dynamics  may  involve  unsteady  and  steady-state  mechanism  which  are 
equilibrium  and  far  from  equilibrium  processes.  In  particular,  studying  the  flow  of  methane  gas,  formation  of 
hydrate,  and  its  dissociation,  along  with  its  complex  mixtures  with  other  hydro-constituents  including 
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sedimentary  compounds  has  attracted  a  considerable  interest  in  marine  geosciences  [9-13].  Of  our  particular 
interest  are  the  effects  of  driving  fields  such  as  concentration  gradients,  hydrostatic  pressure,  thermal  gradient, 
etc.  [14,15]  which  are  crucial  in  understanding  the  hydrate  formation,  mound  evolution,  and  dissociation 
leading  to  such  instabilities  [1]  as  mud  volcanoes  [16-18]  below  the  ocean  floor. 

Enormous  efforts  [19]  have  been  made  in  recent  years  to  understand  mud  volcanoes  [20],  a  violent  eruption 
of  mud  or  clay  often  accompanied  by  methane  gas.  The  main  cause  is  believed  to  be  the  over-pressure  of  the 
fluidized  mud  (clay,  sand,  salt  water)  along  with  a  gas  source  at  depth  within  the  sediment  column  (due  to 
weight  of  mud,  escaping  methane  gas,  etc.).  Most  numerical  modeling  [20]  of  such  complex  phenomena 
involve  balancing  between  the  weight  and  driving  pressure  and  a  source  of  volume  flux  of  simplified  systems. 
Corresponding  competing  terms  offer  good  estimates  of  the  height  reached  by  the  mud  volcanoes  in  the 
framework  of  mean-field  or  first  hand  analysis.  Spatial  heterogeneity  and  local  dynamics  leading  to  such 
cooperative  eruption  seems  to  be  missing  in  such  study.  It  would  be  interesting  to  develop  a  model  where  such 
a  violent  global  response  as  eruption  appears  as  a  collective  cooperative  phenomena  of  its  interacting 
components.  Characterizing  different  components  of  multi-constituents  fluid  mixture  such  as  gas,  liquid,  and 
sediment  by  particles  with  desirable  characteristics  via  coarse  grained  models  seems  a  viable  approach  to 
investigate  such  driven  systems.  In  a  relatively  dense  system  of  such  mixture,  the  dynamics  of  individual 
constituents  and  their  distribution  lead  to  correlations  induced  by  their  interactions.  While  addressing  the 
global  flow  is  the  main  issue  to  address,  the  underlying  self-organizing  structural  morphology  cannot  be 
ignored. 

Investigations  of  structures  and  patterns  in  driven  granular  systems  [21,22]  have  attracted  a  considerable 
interests  in  recent  years  [23-29].  It  is  often  difficult  to  monitor  the  mobility  of  individual  constituents  in 
laboratory  experiments  [30]  in  order  to  understand  their  global  patterns.  On  the  other  hand,  it  is  easier  to 
monitor  their  movements  and  structural  patterns  [31,32]  by  computer  simulations  such  as  lattice  gas  [33,34], 
molecular  dynamics  [26,35,36],  and  Monte  Carlo  [37]  methods.  The  number  of  constituents  is  conserved  in 
most  studies  of  driven  particles  systems  [38,39].  The  self-organizing  system  of  our  driven  mixture  is  non¬ 
conservative  where  the  particles  are  continuously  released  from  a  source  (see  below).  It  should  be  pointed  out 
that  it  is  easier  to  incorporate  interactions  between  constituent  particles  and  effective  medium  (empty  lattice 
sites,  see  below)  in  interacting  lattice  gas  than  that  with  the  Boltzmann  lattice  gas  [40,41]. 

When  an  immiscible  fluid  mixture  consisting  of  two  components  say,  A  and  B  with  dissimilar  masses  (i.e., 
gas  and  liquid,  oil  and  gas,  oil  and  water),  is  released  from  a  source  at  the  bottom  (below  the  ocean  floor), 
several  interesting  questions  arise  regarding  their  flow  and  structure.  How  long  will  it  take  for  constituents  to 
emerge  from  the  top?  How  does  the  morphology  emerge  as  their  sedimentation  occur  at  different  rates?  How 
do  they  self-organize,  mix,  or  segregate?  How  can  structural  evolution  and  the  flow  rate  respond  to  an  external 
hydrostatic  pressure  bias?  As  pointed  out  above,  understanding  such  basic  questions  is  essential  in  context  to 
the  flow  of  methane  gas,  hydrocarbon,  their  complex  compounds,  and  hydrate  formation  and  its  dissociation 
below  the  ocean  floor  [9-13].  Using  an  interacting  lattice  gas  model,  we  have  studied  [42]  the  onset  of 
segregation  and  partial  layering  in  a  steady-state  flow  of  an  immiscible  mixture.  Effects  of  hydrostatic  pressure 
bias  {H)  and  molecular  weight  (a  =  Mb/Ma  =  1-10)  on  steady-state  density  profiles  have  been  recently 
examined  and  attempts  are  made  to  quantify  the  dependence  of  their  overall  density  on  these  parameters 
(//,a)  [43].  Very  recently,  we  studied  the  response  of  flux  density  to  hydrostatic  pressure  bias  [44]  in  a  self¬ 
organizing  immiscible  mixture  where  both  components  (/I,  B)  are  released  into  the  system  with  equal  rate.  We 
observe  linear  response  for  both  constituents  with  rates  depending  on  their  molecular  weights,  especially  at 
low  bias.  We  extend  this  study  here  to  a  system  where  the  release  of  constituents  from  the  source  depends  on 
their  concentration  in  the  system.  Consequently,  we  find  a  dramatic  change  in  response  of  the  flux  density  to 
bias,  from  linear  to  eruptive  (see  below).  The  model  is  described  in  the  next  section  followed  by  results  and 
discussion  with  a  summary  and  conclusion. 

2.  Model  and  method 

We  consider  a  cubic  lattice  of  size  L?  with  L  =  30—200.  Bottom  of  the  lattice  (z  =  1)  is  connected  to  a 
source  releasing  particles  A  and  B  into  the  lattice.  Top  of  the  lattice  (z  =  L)  is  open.  Initially,  about  half  of  the 
lattice  sites  are  randomly  occupied  by  ^4  or  ^  with  equal  probability,  the  remaining  empty  sites  represent  the 
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background  solvent  (S).  The  steady-state  global  properties  do  not  depend  on  the  initial  configuration  (see 
below).  Depending  on  the  occupancy  {A,B,  or  5),  a  lattice  site  has  three  states.  Our  system  with  mobile 
particles  is  characterized  by  an  interaction  energy  E  between  these  components, 

£  =  ^^7(/,A:),  (1) 

/  k 

where  index  /  runs  over  all  particles  and  k  over  all  nearest  neighbor  sites  of  /  with, 

JiA,A)  =  =  -J(A,B)  =  -J{B,A)  =  -£i,  (2) 

J{A,S)  =  J{B,S)=-e2.  (3) 

The  miscibility,  phase  separation  between  A  and  B,  and  their  mobility  are  controlled  by  s\  and  £2.  A  positive 
miscibility  gap  ei  =  1,  considered  here,  enhances  the  phase  separation  in  contrast  to  mixing  with  a  negative  e\. 
The  interaction  with  the  effective  host  medium  is  governed  by  £2  which  is  set  at  ^2  =  1.  Although  these  host 
empty  sites  appear  static,  they  do  exchange  their  positions  with  the  mobile  constituents  (see  below)  and 
therefore,  can  be  used  to  vary  the  viscous  drag  or  the  quality  of  the  solvent  by  varying  the  interaction  62- 

Effect  of  molecular  weights  of  constituents  (M a,  Mb)  is  considered  via  change  in  the  gravitational  potential 
energy  as  they  move  along  (— z)  or  against  (-fz)  gravity.  Further,  a  hydrostatic  pressure  bias  (//),  is  also 
included  to  drive  the  constituents  upward.  The  Metropolis  algorithm  is  used  to  move  randomly  selected 
particles  to  their  neighboring  sites  along  ±x,  ±y,  ±z  directions  with  probabilities, 

P,  =  P_,=Py=  P.y=J-,  P,  =1^,  />_,  =1^;  0^H^\.  (4) 

DO  O 

Change  in  configurational  interaction  energy  along  with  the  bias,  thus,  orchestrate  the  mobility.  Periodic 
boundary  condition  is  used  along  the  transverse  direction.  Top  (z  =  L)  and  bottom  (z  =  1)  of  the  sample  are 
open  for  the  lattice  particles  to  escape.  A  particle  is  released  into  a  bottom  lattice  site  as  soon  as  it  is  vacated  by 
the  particle  when  it  moves  up  or  diffuses  laterally. 

Two  methods  are  considered  to  release  the  type  of  particles  (A  or  B)  from  the  source:  injection  probabilities 
of  A  and  B  are  equal  to  their  relative  concentrations  in  the  lattice  in  procedure  (i)  while  the  probabilities  of 
releasing  particles  (A,  B)  are  equal  in  procedure  (ii)  regardless  of  their  concentrations.  Method  (i)  introduces  a 
stronger  interaction  between  the  source  and  the  sample  which  may  be  more  realistic  in  many  situations,  i.e., 
constituents  of  a  highly  viscous  fluid  have  strong  attractive  interactions  in  general.  Attempt  to  move  each 
particles  once  defines  unit  Monte  Carlo  step  (MCS)  time.  As  the  simulation  proceeds,  the  constituents  move  in 
and  out  of  the  sample.  The  number  of  particles  and  their  density  profiles  (planer  densities  along  x,y,z)  evolve 
but  eventually  stabilize  to  a  steady-state.  Simulations  are  performed  with  many  independent  samples  each  for 
a  sufficiently  long  time  steps  to  reach  steady-state  to  estimate  the  physical  quantities  reliably.  We  have 
observed  [32]  evolution  of  interesting  structural  patterns,  i.e.,  onset  of  phase  separation,  oscillation  in  lateral 
density  profiles,  layering  etc.  which  share  some  characteristics  observed  in  granular  systems  [26].  In  the 
following,  we  focus  mainly  on  the  global  transport  and  flow. 

3.  Results  and  discussion 

Simulations  are  performed  on  samples  of  various  sizes  30*^  to  200^  as  mentioned  above.  However,  data  for 
most  production  runs  are  generated  with  30^,  50^  and  100^  lattices  using  as  many  as  128  independent  runs  with 
largest  sample  for  the  entire  range  of  pressure  bias  H  =  0.0- 1.0.  Finite  size  effects  on  the  qualitative  behavior 
of  the  physical  quantities  are  negligible  particularly  with  50'^  and  100^  samples.  We  restrict  to  miscibility  gap 
£1  =  1.0,  effective  medium  (viscous)  interaction  £2  =  1.0,  molecular  weights  =  0.1,m^  =  0.3,  and 
temperature  T  =  1.0  in  appropriate  unit  of  interaction  energy  and  Boltzmann  constant.  Self-organizing 
structural  evolution  leading  their  steady-state  equilibrium  has  been  already  examined  in  depth  [42,43]  with 
their  different  phases,  segregation,  and  mixing. 

Transport  and  mobility  of  constituents  (tracers)  can  be  analyzed  by  examining  the  variation  of  the  root 
mean  square  (rms)  displacements  of  A  and  B.  Obviously,  the  particles  in  the  dense  phase  (toward  bottom)  are 
not  expected  to  be  as  mobile  as  those  in  the  gas  phase  (toward  top).  Most  particles  entering  at  the  bottom 
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leave  the  sample  from  top  in  the  gas  phase.  The  fast  movements  in  the  gas  phase  dominate  the  asymptotic 
dependence  of  tracers  rms  displacements  on  time  steps.  However,  in  order  to  estimate  the  global  transport,  it  is 
worth  looking  at  the  variation  of  the  rms  displacements  {R)  of  the  center  of  mass  of  each  species  {A  and  B) 
with  the  time  (/)  steps.  Fig.  1  shows  R  versus  t  plots  for  the  various  values  of  bias  H.  The  asymptotic 
dependence  can  be  described  by  a  power-law, 

R  =  Dt\  (5) 


where  D  is  a  constant  and  v,  the  power-law  exponent.  For  most  values  of  the  bias,  the  overall  motion  of  each 
constituent  (A,B)  appears  to  be  drift-like  with  v  ~  1.  At  extremely  low  pressure  bias,  we  see  the  signature  of 
sub-diffusive  motion  (v<0.5  (the  diffusion  movement)). 

Let  us  examine  factors  affecting  the  transport.  The  stochastic  thermal  motion  (the  interplay  between  the 
interaction  energy  and  temperature,  i.e.,  diffusion  in  lateral  (transverse)  planes)  is  augmented  by  three  driving 
fields:  (a)  the  concentration  gradient  caused  by  continuous  release  of  particles  at  the  bottom,  (b)  gravity 
pulling  constituents  downward,  and  (c)  the  external  hydrostatic  pressure  bias  pushing  them  upward.  At  low 
pressure  bias,  the  competition  between  gravity  (downward  movement)  and  the  concentration  gradient 
(upward  from  the  bottom),  appears  to  be  major  cause  for  the  slow  transport  observed  above.  Increasing  the 
bias  offsets  the  gravity,  leading  to  a  net  upward  global  bias  (effect  of  both  concentration  gradient  and 
pressure)  which  drives  the  constituents  into  a  drift-like  motion. 

Because  of  the  net  driving  field  along  the  upward  direction,  there  is  a  net  mass  transfer:  constituents  enter  at 
the  bottom  and  leave  mostly  from  the  top  (see  below).  In  steady-state,  the  net  amount  of  components  entering 
the  system  must  be  equal  to  net  mass  leaving  the  system,  i.e.,  the  conservation  of  mass  leading  to  the 
continuity  equation.  From  the  variation  of  the  net  mass  transfer  (Q)  in  time  /,  we  can  estimate  the  current  due 
to  mass  flow. 


(6) 


The  current  density  (/),  current  per  unit  cross-sectional  area  {Ac), 


J  =  i/Ac- 


(7) 
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Fig.  1.  RMS  displacements  of  the  center  of  mass  of  A  and  B  versus  time  t  steps  at  various  values  of  H  with  f.\  =  \.  Sample  size  100^  with 
64-128  independent  runs  are  used. 


420 


R.B.  Pandey,  J.F.  Get  trust  /  Physica  A  368  (2006)  416-424 


In  presence  of  a  field  A/?  (pressure,  flow  field,  concentration,  etc.  and  their  combination),  one  may  expect,  at 
least  in  simple  cases,  a  power-law  dependence, 

j  =  (8) 

where  a  is  an  analog  of  conductivity.  The  exponent  (i  =  1  corresponds  to  the  well-known  Darcy’s  law  of 
linear  response.  In  our  case,  it  is  rather  difficult  to  estimate  Ap  due  to  different  driving  mechanisms  mentioned 
above.  However,  we  can  examine  our  data  to  extract  the  response  of  j  to  pressure  bias  (//),  the  variable  in  our 
study. 

According  to  conservation  of  mass  (continuity  equation),  the  flux  current  density  (/)  must  be  independent  of 
time  in  the  asymptotic  regime.  Fig.  2  shows  the  temporal  variation  of  j  for  B  (heavier)  constituents  (similar 
asymptotic  behavior  is  also  observed  for  A  constituents).  We  see  that  the  flux  densities  at  both  bottom  and  top 
become  equal  after  an  initial  relaxation  time  and  remain  constant  thereafter.  The  current  density  depends  on 
the  pressure  bias  (//).  Variations  of  and  jg  of  constituents  A  and  B,  respectively,  with  H  are  presented  in 
Figs.  3  and  4.  A  linear  response  is  generally  seen  in  both  systems,  i.e.,  with  concentration  dependent  (i)  (Fig.  3) 
and  independent  (ii)  (Fig.  4),  for  most  region  of  H  which  implies  the  validity  of  the  Darcy’s  law.  When 
constituents  A  and  B  are  released  with  equal  probabilities  (ii),  their  respective  current  densities  respond 

linearly  to  nearly  entire  range  of  bias  H  (i.e.,  i?  ~  1);  the  prefactor  a  is  larger  for  the  heavier  component  (B). 
The  response  of  current  density  for  system  (i),  on  the  other  hand,  is  much  more  complex.  The  flow  of  A 
responds  differently  than  that  of  B  as  expected  due  to  difference  in  their  molecular  weights.  While  the  flow 
density  of  constituents  A  is  larger  than  jg  at  lower  bias  (//^0.4),  the  flow  of  constituents  B  becomes  larger 
than  that  of  A,  i.e.,y^  <jg,  at  higher  bias  (//  ^0.4).  Since  B  has  larger  molecular  weight,  they  tend  to  gravitate 
more  than  A  toward  the  bottom.  As  a  result,  particles  B  have  higher  probability  of  leaking  from  the  bottom 
than  particles  A  particularly  at  lower  bias.  Increasing  the  upward  bias  (//)  reduces  the  effect  of  gravitation 
which  results  in  a  larger  jg.  Relatively  larger  amount  of  B  constituents  enter  the  sample  at  larger  H.  Thus,  the 
competition  between  gravity  and  hydrostatic  pressure  bias,  leads  to  the  crossover  from  A  dominated  flow 
Ua  ^Jb)  lower  bias  H  to  B  dominated  flow  (J^  <jg)  at  higher  H. 

At  extreme  values  of  bias,  the  response  becomes  non-linear  with  a  rapid  increasing  response  of  B  leading  to 
a  “eruptive  flow’’.  At  the  same  time,  the  flow  response  of  A  becomes  non-monotonic  with  a  rapid  decay  of  j^ 


Fig.  2.  Flux  density  j  of  B  from  top  (open)  and  bottom  (filled)  with  time  steps  for  the  same  parameters  as  in  Fig.  1. 
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Both  systems  (i)  and  (ii)  are  open  for  constituents  to  escape.  There  is  a  net  flow  of  mass  of  these  immiscible 
constituents  (A,B)  without  conserving  their  concentrations.  However,  steady-state  density  profiles  are 
obtained  and  the  continuity  equations  for  the  conservation  of  mass  of  each  components  are  satisfied  within  the 
range  of  statistical  fluctuations.  If  the  release  of  particles  from  the  source  is  independent  of  their 
concentration,  we  find  that  corresponding  current  densities  respond  linearly  to  bias  (//).  Further, 

Ja^Jb  except  at  extreme  value  of  H  (->  1),  —Jb’  system  /  where  the  reservoir  and  the  lattice  arc 
correlated  with  a  concentration  dependent  release  rate  of  constituents,  the  flow  response  is  complex.  We 
observe  a  crossover  in  response  of  the  current  densities  from  low  region  of  bias  (//^0.4)  with  >Jb  to  higher 
range  of  bias  (//  >0.4)  withy^  <y^.  Constituents  (A)  with  lower  molecular  weight  respond  linearly  to  bias  for 
most  region  of  H  except  at  high  values  (//^0.8)  where  the  response  becomes  non-monotonic.  The  heavier 
component  {B)  responds  nonlinearly  on  increasing  the  bias:  a  large  response  at  low  values  of  H  followed  by 
linear  response  before  an  vulcanizing  response  in  high  range  of  H.  We  have  defined  a  volatility  response 
parameter  Vb  which  seems  to  diverge  ->  oo  as  //  ->  1;  corresponding  parameter  for  A  becomes  negative 
in  this  range. 

In  summary,  we  have  identified  the  linear  response  region  of  both  immiscible  components  (A,B)  where 
Darcy's  law  is  valid  and  provided  empirical  relations.  Differences  in  response  of  two  immiscible  components 
depending  on  their  molecular  weight  are  identified  in  two  systems  of  source  reservoirs.  We  are  able  to  identify 
the  condition  in  which  an  eruptive  flow  response  occurs  which  may  complement  existing  modeling  [20]  and 
stimulate  further  interest  in  this  field  [19]. 
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